########################################################################################################################
## This program plots the graph of estimates and counterfactuals ##
## Model refers to "Distribution Regression with Selection"
## Authors: Chernozhukov, V., I. Fernandez-Val and S. Luo
## Last Update: 6/17/2018

## 80s vs 90s vs 100s within gender 

## male vs female within year(s)




rm(list=ls());

library(fields);


#####################
# Extract Arguments #
#####################
#Define default values for variables
counter.method <- 2; 
# 1 for counterfactual decades within gender; 
# 2 for counterfactual sex within time period; 
# 3 for counterfactual 18 years within gender;
rho.spec       <- 1; 
# Specification of rho: 1 for default; 
#                       2 for education-specific rho's; 
#                       3 for couple-specific rho's;
#                       4 for time-specific rho's;
#                       5 for linear expression of X;
#                       6 for age; 
#                       7 for number of kids;
#                       8 for linear expression of years of schooling;
#                       9 for 4 education dummies;
#                       10 for 4 education dummies x couple;
#                       11 for age x couple;
#                       12 for 4 age dummies;
#                       13 for numkids x couple;
#                       14 for 4 age dummies x couple;
#                       15 for linear in t;
#                       16 for quardratic in t;
#                       17 for yrseduc x couple;
#                       18 for linear in t x couple;
#                       19 for 6 time dummies;
#                       20 for propensity score of employment;
#                       21 for pscore x couple;
sample.id      <- 11378;
l.trim         <- 10;
u.trim         <- 90;
alpha          <- 0.05;
n.boot         <- 500;


str.sample <- ifelse(counter.method==1, ifelse(sample.id==1, "_m","_f"), 
                     ifelse(counter.method==2, paste("_",sample.id,sep=""), ifelse(sample.id==1, "_mhalf","_fhalf")));
str.rho.all <- c("","_educ","_couple","_time","_allx","_age","_numkids","_yrseduc","_4dum","_4dumxcouple",
                 "_agexcouple","_4agedum","_numkidsxcouple","_4agedumxcouple","_lintime","_quartime","_yrseducxcouple",
                 "_lintimexcouple","_6tdum","_pscore","_pscorexcouple");
str.rho <- str.rho.all[rho.spec];

#str.trim <- paste("_", ifelse(l.trim<10, "0",""), l.trim, ifelse(u.trim<10, "0",""), u.trim, sep="");
str.trim <- paste("_", l.trim, u.trim, sep="")

str.cb <- ifelse(alpha==0.05, "","_c90");

str.robust <- "";

str.boot <- ifelse(n.boot==200,"",paste("_",n.boot,sep=""));


##########################
# Load Inference results #
##########################
load(paste("Results",str.sample,str.rho,str.trim,str.cb,str.robust,str.boot,".RData",sep=""));


qn.plot <- c(qn.lower:qn.upper);
q.plot <- seq((l.thres+(qn.lower-1)*by.thres),(l.thres+(qn.upper-1)*by.thres),by.thres)*100;
if(counter.method == 1){
  main.str <- ifelse(sample.id==1, "Male","Female");
  main.str.counter1 <- paste(main.str," in 1980's", sep="");
  main.str.counter2 <- paste(main.str," in 2000's", sep="");
  
  legend.counter1 <- "1980s";
  legend.counter2 <- "2000s";
}else if(counter.method == 3){
  main.str <- ifelse(sample.id==1, " - Male"," - Female");
  main.str.counter1 <- paste(main.str," in recent half: 1996-2013", sep="");
  main.str.counter2 <- paste(main.str," in previous half: 1978-1995", sep="");
  
  legend.counter1 <- "recent half";
  legend.counter2 <- "prev half";
}else{
  if (sample.id >= 1000){
    if(sample.id >= 100000){
      from.year <- sample.id %/% 1000;
      to.year   <- sample.id - from.year*1000;
    }else{
      from.year <- sample.id %/% 100;
      to.year   <- sample.id - from.year*100;
    };
    to.year.str <- 1900 + to.year;
    from.year.str <- 1900 + from.year;
    
    main.str <- "";
    #main.str.counter1 <- paste("Male in ", to.year.str, " ~ ", from.year.str, sep="");
    #main.str.counter2 <- paste("Female in ", to.year.str, " ~ ", from.year.str, sep="");
    main.str.counter1 <- "Male";
    main.str.counter2 <- "Female";
  }else{
    sample.id.str <- 1900 + sample.id;
    main.str <- sample.id.str;
    main.str.counter1 <- paste("Male in ", sample.id.str, sep="");
    main.str.counter2 <- paste("Female in ", sample.id.str, sep="");
  };
  
  legend.counter1 <- "male";
  legend.counter2 <- "female";
};



############################
# First Stage Result Table #
############################
fstable <- cbind(c(rbind(colnames(elist.counter1$pi)[c(1:11,23:28,(NCOL(elist.counter1$pi)-1):NCOL(elist.counter1$pi))],rep("",19))),
                 c(rbind(format(elist.counter1$pi[1,c(1:11,23:28,(NCOL(elist.counter1$pi)-1):NCOL(elist.counter1$pi))],digits=3),
                         paste("(",format(elist.counter1$pi[2,c(1:11,23:28,(NCOL(elist.counter1$pi)-1):NCOL(elist.counter1$pi))],digits=2),")",sep=""))),
                 c(rbind(format(elist.counter2$pi[1,c(1:11,23:28,(NCOL(elist.counter1$pi)-1):NCOL(elist.counter1$pi))],digits=3),
                         paste("(",format(elist.counter2$pi[2,c(1:11,23:28,(NCOL(elist.counter1$pi)-1):NCOL(elist.counter1$pi))],digits=2),")",sep=""))));
colnames(fstable) <- c("Variables","Male","Female");

library(xtable);
options(xtable.floating=FALSE);
options(xtable.timestamp="");

print(xtable(fstable),include.rownames=F);


################
# Plot Setting #
################
para.namelist <- c("educ16","educ1718","educ1920","educ2122","educ23","couple");
#c("educ16","educ1718","educ1920","educ2122","educ23","couple"); # "educ23",;
# "agep","agep2","agep3","agep4",
# "factor(region)2","factor(region)3","factor(region)4","factor(region)5","factor(region)6",
# "factor(region)7","factor(region)8","factor(region)9","factor(region)10","factor(region)11",
# "factor(region)12","numch1","numch2","numch34","numch510","numch1116","numch1718");
ylimlist      <- rbind(c(0,2.5),c(0,2.5),c(0,2.5),c(0,2.5),c(0,2.5),c(-0.3,0.7));

if(rho.spec == 2){
  delta.each <- c("educ15","educ16","educ1718","educ1920","educ2122","educ23");
  mfrow.current <- c(3,2);
}else if(rho.spec == 3){
  delta.each <- c("single","couple");
  mfrow.current <- c(2,1);
}else if(rho.spec == 4 | rho.spec == 19){
  delta.each <- colnames(model.matrix(rho.formula, mysample));
  mfrow.current <- c(3,2);
}else if(rho.spec == 9){
  delta.each <- c("hsdrop","hs","somecol","col");
  mfrow.current <- c(2,2);
}else if(rho.spec == 10){
  delta.each <- colnames(model.matrix(rho.formula, mysample));
  mfrow.current <- c(3,3);
}else if(rho.spec == 12){
  delta.each <- colnames(model.matrix(rho.formula, mysample));
  mfrow.current <- c(2,2);
}else if(rho.spec == 14){
  delta.each <- colnames(model.matrix(rho.formula, mysample));
  mfrow.current <- c(3,3);
};

if(rho.spec == 5){
  delta.namelist <- c("educ16","educ1718","educ1920","educ2122","educ23","couple");
  mfrow.current  <- c(3,2);
  ylim.current <- c(-3,4);
}else if(rho.spec == 6){
  delta.namelist <- "age";
  mfrow.current  <- c(1,1);
  ylim.current <- c(-0.1,0.1);
}else if(rho.spec == 7){
  delta.namelist <- c("numch1","numch2","numch34","numch510","numch1116","numch1718");
  mfrow.current  <- c(3,2);
  ylim.current <- c(-3,4);
}else if(rho.spec == 8){
  delta.namelist <- "yrseduc";
  mfrow.current  <- c(1,1);
  ylim.current <- c(-0.1,0.1);
}else if(rho.spec == 11){
  delta.namelist <- c("age","couple:age");
  mfrow.current <- c(2,1);
  ylim.current <- c(-0.1,0.1);
}else if(rho.spec == 13){
  delta.namelist <- c("numchyt6","couple:numchyt6");
  mfrow.current <- c(2,1);
  ylim.current <- c(-0.7,1);
}else if(rho.spec == 15){
  delta.namelist <- c("(Intercept)","year_res");
  mfrow.current <- c(2,1);
}else if(rho.spec == 16){
  delta.namelist <- c("(Intercept)","year_res","I(year_res^2)");
  mfrow.current <- c(2,2);
}else if(rho.spec == 17){
  delta.namelist <- c("(Intercept)","yrseduc","couple","yrseduc:couple");
  mfrow.current <- c(2,2);
  ylim.current <- c(-1,1);
}else if(rho.spec == 18){
  delta.namelist <- c("(Intercept)","couple","year_res","year_res:couple");
  mfrow.current <- c(2,2);
}else if(rho.spec == 20){
  delta.namelist <- c("(Intercept)","pscore");
  mfrow.current <- c(2,1);
  ylim.current <- c(-2,2);
}else if(rho.spec == 21){
  delta.namelist <- c("(Intercept)","pscore","couple","pscore:couple");
  mfrow.current <- c(2,2);
  ylim.current <- c(-3,3);
};


if(counter.method==2){
  ylim.pct.F <- c(-50,130);
  ylim.pct.Fcond <- c(-50,130);
  ylim.QF.diff <- c(-0.2,0.5);
  ylim.QF.cond.diff <- c(-0.3,0.5);
}else if(counter.method == 1 | counter.method==3){
  if(sample.id==1){
    ylim.pct.F <- c(-100,150);
    ylim.pct.Fcond <- c(-50,150);
    ylim.QF.diff <- c(-0.4,0.6);
    ylim.QF.cond.diff <- c(-0.4,0.6);
  }else{
    ylim.pct.F <- c(-50,150);
    ylim.pct.Fcond <- c(-50,150);
    ylim.QF.diff <- c(-0.2,0.6);
    ylim.QF.cond.diff <- c(-0.2,0.5);
  }
};



source("DRS_Emp_Plot_Para.R");

source("DRS_Emp_Plot_DFQF.R");

source("DRS_Emp_Plot_Cond.R");
